I'm gonna port this to a script soon. I need to constrain the RedMagic HOD against nd, fc, and $\langle N{gal} | M_{>e14} \rangle$. I don't need to do any populations for this so it should be quick.


In [49]:
from pearce.mocks.kittens import cat_dict
import numpy as np
from astropy.cosmology import LambdaCDM
from astropy.io import fits
from scipy.linalg import inv
from os import path

In [50]:
a = 0.81120
z = 1./a -1
cosmo_params = {'simname':'chinchilla', 'Lbox':400.0, 'scale_factors':[a]}

In [51]:
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!

In [52]:
cat.load(a, tol = 0.01, HOD='redMagic', particles = False)#, hod_kwargs = {'sec_haloprop_key':'halo_log_nfw_conc'})#, hod_kwargs={'split': 0.5})

In [53]:
#vals to constrain

fname = '/u/ki/jderose/public_html/bcc/measurement/y3/3x2pt/buzzard/flock/buzzard-2/tpt_Y3_v0.fits'
hdulist = fits.open(fname)
zbin = 1
z_bins = np.array([0.15, 0.3, 0.45, 0.6, 0.75, 0.9])
nz_zspec = hdulist[8]
N= sum(row[2+zbin] for row in nz_zspec.data)

area = 5063 #sq degrees
full_sky = 41253 #sq degrees

buzzard = LambdaCDM(H0=70, Om0=0.286, Ode0=0.714, Tcmb0=2.725, Neff=3.04)
#volIn, volOut = buzzard.comoving_volume(z_bins[zbin-1]), buzzard.comoving_volume(z_bins[zbin])
volIn, volOut = buzzard.comoving_volume(z_bins[zbin-1]), buzzard.comoving_volume(z_bins[zbin])

fullsky_volume = volOut-volIn
survey_volume = fullsky_volume*area/full_sky
nd = N/survey_volume
print nd
nd_std = np.sqrt(N)/survey_volume
print nd_std


0.000284539736705 1 / Mpc3
6.11545845728e-07 1 / Mpc3

In [54]:
fc = 0.2
fc_std = 0.1

In [55]:
Nm14 = 4
Nmstd = 2

In [56]:
param_names = cat.model.param_dict.keys()
param_bounds = {'logMmin': [11.5, 13.5], 'sigma_logM': [0.05, 0.6], 'logM0': [12.0, 16.0],\
                'logM1': [13.0, 15.0], 'alpha': [0.8, 1.2], 'f_c': [0.01, 0.5]}

In [67]:
pname = 'sigma_logM'
vals = np.linspace(param_bounds[pname][0], param_bounds[pname][1], 10)

params = {'logMmin': 12.1, 'sigma_logM': 0.3, 'logM0': 14.0,\
                'logM1': 14.0, 'alpha': 1.0, 'f_c': 0.19}

for v in vals:
    params[pname] = v
    print v, cat.calc_analytic_nd(params)


 0.05 0.00020628909633046424
0.1111111111111111 0.00020851049937752077
0.17222222222222222 0.00021247207162946357
0.23333333333333328 0.00021823956539403035
0.2944444444444444 0.00022592997789618545
0.3555555555555555 0.00023571819784037686
0.4166666666666666 0.00024778132252039324
0.4777777777777777 0.00026215715888630017
0.5388888888888889 0.00027861799966679974
0.6 0.00029669791593415535

In [57]:
cov = np.diag(np.array([nd_std.value**2, fc_std**2, Nmstd**2]) )
invcov = inv(cov)

In [58]:
def lnprior(theta, param_names, param_bounds, *args):
    for p, t in izip(param_names, theta):
        low, high = param_bounds[p]

        if np.isnan(t) or t < low or t > high:
            return -np.inf
    return 0

In [59]:
hod_kwargs = {'mass_bin_range': (9,16),
'mass_bin_size': 0.01,
'min_ptcl': 200}
mf = cat.calc_mf(**hod_kwargs)
mass_bins = np.logspace(hod_kwargs['mass_bin_range'][0],\
                        hod_kwargs['mass_bin_range'][1],\
                        int( (hod_kwargs['mass_bin_range'][1]-hod_kwargs['mass_bin_range'][0])/hod_kwargs['mass_bin_size'] )+1 )

In [60]:
def lnlike(theta, param_names, param_bounds, obs_vals, invcov, mf, mass_bins, hod_kwargs):
    params = dict(zip(param_names, theta))
    f_c = params['f_c']
    
    hod = self.calc_hod(params, **hod_kwargs)
    nd = np.sum(mf*hod)/((self.Lbox/self.h)**3)
    Nm14 = np.mean(hod[mass_bins>10**14])
    
    pred_vals = np.array([nd, f_c, Nm14])
    delta = pred_vals - obs_vals
    return -delta.dot(invcov.dot(delta))

In [61]:
def lnprob(theta, *args):
    """
    The total liklihood for an MCMC. Mostly a generic wrapper for the below functions.
    :param theta:
        Parameters for the proposal
    :param args:
        Arguments to pass into the liklihood
    :return:
        Log Liklihood of theta, a float.
    """
    lp = lnprior(theta, *args)
    if not np.isfinite(lp):
        return -np.inf

    return lp + lnlike(theta, *args)

In [62]:
nwalkers = 500
nsteps = 5000
nburn = 0

In [63]:
savedir = '/home/users/swmclau2/scratch/PearceMCMC/'
chain_fname = path.join(savedir,'%d_walkers_%d_steps_chain_wt_alt_redmagic_z%.2f.npy'%(nwalkers, nsteps, z))

with open(chain_fname, 'w') as f:
    f.write('#' + '\t'.join(param_names)+'\n')


---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-63-b21bfef08370> in <module>()
      2 chain_fname = path.join(savedir,'%d_walkers_%d_steps_chain_wt_alt_redmagic_z%.2f.npy'%(nwalkers, nsteps, z))
      3 
----> 4 with open(chain_fname, 'w') as f:
      5     f.write('#' + '\t'.join(param_names)+'\n')

IOError: [Errno 2] No such file or directory: '/home/users/swmclau2/scratch/PearceMCMC/500_walkers_5000_steps_chain_wt_alt_redmagic_z0.23.npy'

In [ ]:
ncores = 1
num_params = len(param_names)

sampler = mc.EnsembleSampler(nwalkers, num_params, lnprob,
                             threads=ncores, args=(param_names, param_bounds, obs_vals, invcov, mf, mass_bins, hod_kwargs))


for result in sampler.sample(pos0, iterations=nsteps, storechain=False):
    with open(chain_fname, 'a') as f:
        np.savetxt(f, pos[0])